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It is well known that the two-dimensional (2D) nonlinear Schrodinger equation (NLSE) with 
the cubic-quintic (CQ) nonlinearity supports a family of stable fundamental solitons, as well as 
solitary vortices (alias vortex rings), which are stable for sufficiently large values of the norm. We 
study stationary localized modes in a symmetric linearly coupled system of two such equations, 
focusing on asymmetric states. The model may describe "optical bullets" in dual-core nonlinear 
optical waveguides (including spatiotemporal vortices that were not discussed before), or a Bose- 
Einstein condensate (BEC) loaded into a "dual-pancake" trap. Each family of solutions in the 
single-component model has two different counterparts in the coupled system, one symmetric and 
one asymmetric. Similarly to the earlier studied coupled ID system with the CQ nonlinearity, 
the present model features bifurcation loops, for fundamental and vortex solitons alike: with the 
increase of the total energy (norm), the symmetric solitons become unstable at a point of the direct 
bifurcation, which is followed, at larger values of the energy, by the reverse bifurcation restabilizing 
the symmetric solitons. However, on the contrary to the ID system, both the direct and reverse 
bifurcation may be of the subcritical type, at sufficiently small values of the coupling constant, A. 
Thus, the system demonstrates a double bistability for the fundamental solitons. The stability of 
the solitons is investigated via the computation of instability growth rates for small perturbations. 
Vortex rings, which we study for two values of the "spin" , s = 1 and 2, may be subject to the 
azimuthal instability, like in the single-component model. In particular, complete destabilization of 
asymmetric vortices is demonstrated for a sufficiently strong linear coupling. With the decrease of 
A, a region of stable asymmetric vortices appears, and a single region of bistability for the vortices is 
found. We also develop a quasi-analytical approach to the description of the bifurcations diagrams, 
based on the variational approximation. Splitting of asymmetric vortices, induced by the azimuthal 
instability, is studied by means of direct simulations. Interactions between initially quiescent solitons 
of different types are studied too. In particular, we confirm the prediction of the reversal of the sign 
of the interaction (attractive/repulsive for in-phase/out-of-phase pairs) for the solitons with the odd 
spin, s = 1, in comparison with the even values, s = and 2. 

PACS numbers: 42.65.Tg; 03.75.Lm; 05.45.Yv; 47.20.Ky 
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I. INTRODUCTION 



In the last two decades, much interest has been drawn to the studies of spatiotemporal sohtons (STSs, ahas "hght 
buhets") in nonhnear optics [HE]. They are supported through the balance between the temporal dispersion, spatial 
diffraction, and nonlinearity. It is well known that two-dimensional (2D) and three-dimensional (3D) STSs cannot be 
stable in uniform media with the cubic (Kerr) nonlinearity, due to possibility of the collapse in the same setting [2l[3]. 
To avoid the collapse, other types of the nonlinearity were proposed. In particular, the stability of multidimensional 
solitons is readily secured by saturable [HIS], quadratic (x*^^^) [SI (Zj: and cubic-quintic (CQ) nonlinearities [8l[9]. 
While the creation of STSs in 3D has not been reported so far, quasi-2D STSs were made in x*^^^ crystals [10^. 

While the fundamental (zero-vorticity) STSs supported by the above-mentioned non-Kerr nonlinearities are stable, 
the stability of solitary vortices (alias vortex rings, or spinning solitons, with integer "spin" s referring to the corre- 
sponding topological charge) are vulnerable to the destabilization by azimuthal perturbations. In particular, spinning 
solitons in 2D models with x^^^ or saturable nonlinearity are unstable, as demonstrated by simulations [11 and the 
experiment [12 . The azimuthal instability breaks the solitary vortices with 5 = 1 into two or three fragments, each 
re-trapping into a moving fundamental soliton, so that the intrinsic spin moment of the unstable mode is transformed 
into the orbital momentum of the set of separating fragments. 

Vortex solitons may be stable in media with competing nonlinearities. In particular, stable 2D vortices, with 5 = 1 
and 2, were found in a model combining the x^^^ and self-defocusing cubic nonlinearities ^j. Another model which 
is known to support stable 2D spinning solitons includes self- focusing cubic and self-defocusing quintic nonlinearities. 
Quiroga-Teixeiro and Michinel [14] were the first to demonstrate in direct simulations that 2D solitons with 5 = 1 
may be stable in the CQ model, provided that their power (norm) is sufficiently large. The detailed analysis [15^,16], 
which made use of the linearized version of the model for computing growth rates of perturbation eigenmodes, had 
demonstrated that, for spins 5 = 1 and 2, the CQ model exhibits relatively large stability regions, which cover, 
respectively, ^ 9% and ^ 8% of the respective existence regions, in terms of soliton's norm (total power). It has also 
been demonstrated that the spinning solitons with the vorticity up to 5 = 5 may also be stable, but in extremely 
narrow regions [TT] . 

The subject of the present work are 2D solitons, both the fundamental and spinning ones, in a symmetric system of 
two linearly coupled nonlinear Schrodinger equations (NLSEs) with the CQ nonlinear terms. The main objectives of 
the analysis are the symmetry-breaking bifurcations and asymmetric solitons generated by them. As explained below, 
the model may have realizations in both nonlinear optics and BECs (Bose-Einstein condensates, where the NLSE is 
known as the Gross-Pitaevskii equation, GPE [E]). The ID version of the model was studied in Ref. [19], where it 
was found that a bifurcation of the supercritical type destabilizes symmetric two-component solitons, generating a 
pair of asymmetric ones. At larger values of the norm (power), the branches of symmetric and asymmetric solitons 
merge back, restabilizing the symmetric ones, which gives rise to a bifurcation loop. At relatively small values of the 
linear-coupling constant. A, the reverse bifurcation is of the subcritical type, which may feature a large region of the 
bistability between the symmetric and asymmetric ID solitons. With the increase of A, the bistability region gradually 
disappears and the bifurcation loop shrinks, vanishing and leaving only the symmetric solitons in the system at still 
larger values of A. 

Spontaneous-symmetry-breaking bifurcations of 2D solitons and vortices in linearly-coupled systems were studied in 
Ref. [20 , in terms of two parallel pancake-shaped BECs linked by the tunneling of atoms across the potential barrier 
separating them. The model was based on two GPEs with the linear coupling and cubic nonlinearity. The stability 
of the solitons and vortices against the collapse and (as concerns the vortices) against the splitting was provided by a 
2D periodic potential (an optical lattice) present in both equations. In agreement with the results known from other 
models, the symmetric solitons and vortices underwent the symmetry-breaking bifurcations in the system with the 
self-attractive nonlinearity. Unlike that setting, in the present work we consider the uniform space, the stability of 
the single-component solitons and vortices being provided by the CQ nonlinearity. 

The paper is organized as follows. The model is introduced, and its physical realizations are discussed, in section [llj 
The set of stationary solutions to the coupled equations is parameterized by the linear-coupling constant. A, and 



the propagation coefficient, k. In section |III| we present asymmetric solutions, and produce bifurcation diagrams for 
the fundamental and spinning solitons, with vorticities s = 0,1,2. Section [TVj reports a detailed stability analysis, 
performed by dint of numerical calculations of the growth rates for eigenmodes of small azimuthal perturbations. In 
section [V| we develop an explicit analytical approximation for the solutions, based on the variational method. Results 
of direct numerical simulations, which demonstrate the splitting of azimuthally unstable asymmetric spinning solitons, 
are displayed in section 



concluded by section V 



VI In section VII we report numerical results for interactions between solitons. The paper is 
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II. THE MODEL 



The system of 2D linearly coupled equations with the CQ nonlinearity is taken in the scaled form (cf. the ID 
system introduced in Ref. [H]): 



^yy 



In the terms of the BEC, this is the system of GPEs for wave functions of the condensate in parallel tunnel-coupled 
pancake- shaped traps, with the negative scattering length accounting for the cubic self-attraction, the scattering 
length itself being eliminated by the rescaling [20 . In this context, evolution variable z is actually time, the quintic 
terms account for repulsive three-body collisions, provided that collision-induced losses may be neglected [21 . 

Actually, Eqs. ([T]) may find a more relevant interpretation in the application to optics, where the equations appear 
as normalized NLSEs for the transmission of spatiotemporal light signals in a dual-core planar waveguides. In this 
context, z is the propagation distance, x is the transverse coordinate, and y is the temporal variable (reduced time), 
provided that the sign of the group- velocity dispersion in the waveguide is anomalous [1 . Accordingly, terms ipxx^ (j^xx 
and ^l^yy^ 4>yy accouut for the paraxial diffraction and dispersion of light, respectively, while 1/A determines the coupling 
length in the dual-core waveguide. The equations are made symmetric with respect to x and y by means of rescaling 
of the spatial and temporal variables. Also, A is fixed to be real and positive, which can also be achieved in the general 
case by means of an obvious transformation. As for the self-focusing-defocusing CQ nonlinearity, it was theoretically 
predicted [22] and observed [23] in diverse optical media. Some of them admit the fabrication of dual-core planar 
waveguides. 

The fundamental STS in a planar waveguide, although it has never been reported in an experiment, is a well-known 
concept [2 . On the other hand, the spatiotemporal vortex^ whose snapshot would seem as an elliptic ring running 
at the speed of light in the plane of the waveguide, is a novel object. In the experiment, it may be coupled into 
the planar waveguide by an oblique vortical laser beam shone onto the waveguide under an appropriate angle. In 
that sense, the physical purport of the spatiotemporal vortices is essentially different from that of (2 + l)-dimensional 
spatial solitons with the embedded vorticity, which are understood as hollow cylindrical beams of light propagating 
in a bulk medium [TTTfTB] . In the experiment, vortical spatial solitons, built as multi-beam complexes with the phase 
distribution carrying the effective vorticity (rather than cylindrical beams), were created in photorefractive crystals 
with the saturable nonlinearity, their stability against splitting being maintained by a photoinduced lattice potential 

EH). 

We aim to find stationary axisymmetric solutions to Eqs. ([T]) as 

iIj = U{r) ex.-p{isO) exp(i/cz), 

(j) = V{r) ex.-p{isO) exp(i/cz), (2) 

where r and 6 are the polar coordinates in the (x^y) plane, k is the propagation constant, and integer s is the 
above-mentioned spin. Substituting expressions ^ into Eqs. ([T]), we arrive at equations for real functions U and V: 

(PU IdU . . 

-kV + + ~ - '^V + -V^ = -\U, (3) 
ar^ r ar 

with the boundary conditions demanding that the solution must feature asymptotic forms r'^' at r 0, and 
exp(— V^r) at r ^ oo (hence k must be positive). The energies (norms) of the two components of the soliton 
(alias their norms) are defined as usual, 

/ + 00 
r{U^V^)dr, (4) 
-CO 

the total energy being £^totai = + ^v- The asymmetry of the two-component soliton is characterized by ratio 

e = f^. (5) 
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FIG. 1. Examples of asymmetric solitons with vorticities s = (a), s = 1 (b), and s = 2 (c), for coupling constant A = 0.05. 
Both components are shown, with values of propagation constant k indicated near the corresponding curves. Notice that in the 
{k, B) plane (G is the asymmetry measure defined as per Eq. ([5|), families of the solutions form closed loops, which include the 
reverse bifurcation of the subcritical type, see Figs. |2|4| below. Therefore, asymmetric solutions are found for k increasing up 
to a certain maximum value, and then turning back and decreasing until hitting the reverse-bifurcation point. In the present 
panels, the solitons pertaining to the "backward" subfamily are labeled by (2), to distinguish them from their counterparts 
with the same values of k belonging to the "forward" part of the family. 



III. ASYMMETRIC SOLITONS AND BIFURCATION LOOPS 



Stationary symmetric and asymmetric soliton solutions for s = 0, 1 and 2 were generated in a numerical form, 
applying the Newton- Raphson method to Eqs. (|3|. We have also examined the possibility of the existence of elliptic 
solitons (i.e., anisotropic localized solutions to Eqs. (|3|), using several numerical algorithms adjusted for the 2D setting, 
such as a generalized Petviashvili iteration method, and a modification of the squared-operator method presented in 
Ref. [25]. No stationary elliptic solutions have been found. 

Symmetric soliton solutions in the present model, with tjj = and = 0, can be obtained from their counterparts 
previously found in the single-component CQ model, i.e., one equation from system ([T]), with A = 0, for a single wave 
function, (p: ilj{x^y, z;k) = (j){x,y^z;k) = (p{x,y^z;k — A). Accordingly, the symmetric solutions emerge at /c = A, 
with the known minimum (threshold) values of the energy in the single component [16]: ^thr = 11.73, 48.38, and 
88.34, for s = 0, 1, and 2, respectively. Within the family of the symmetric solitons, k varies from kmin = A, which 
corresponds to E = £^thr, up to /Cmax = A + 3/16, corresponding to £^ ^ oo (/c = 3/16 is the value of the propagation 
constant in the single-component 2D model at which the energy diverges, along with the soliton's radius, for any s 

m)- 

For each value of the spin, a unique family of asymmetric solitons, with U ^ V and O ^ 0, can be found. Typical 
radial profiles of asymmetric solitons with s = 0, 1 and 2, for A = 0.05 and several values of /c, are shown in Fig. [ij 

Similar to its ID counterpart [l^, the present system features bifurcation loops accounting for the transition from 
symmetric solitons to the asymmetric ones and back. Several generic examples of the bifurcation loops are shown in 



5 



Figs. [2j [3]and[4j for 5 = 0, 1 and 2, respectively. As in the ID case, the loops shrink as the coupling constant, A, 
increases, and they expand as A decreases. In the limit of A ^ 0, the loops open up in the direction of £^ ^ 00. 

For 5 = 0, the loop collapses and disappears at Amax^^ ^ 0.0964. Up to A ~ 0.0852, both the direct bifurcation 
and the reverse one, which closes the loop, are subcritical, giving rise to two regions of bistability (which may also 
be called tristability, as the asymmetric soliton always exists in two copies, which are specular images to each other). 
This bifurcation picture is different from the one obtained in the ID model, in which the direct bifurcation is always 
supercritical. In the narrow interval of 0.0852 < A < 0.0891, the direct bifurcation in the present model is supercritical, 
while the reverse one remains subcritical. With the further increase of A up to the point of the disappearance of the 
loop, 0.0891 < A < 0.0964, both the direct and reverse bifurcations are supercritical, and the loop's shape is completely 
convex, featuring no bistability. 

The stability of all the branches of the fundamental {s = 0) solitons strictly follows criteria of the elementary 
bifurcation theory [27]. In particular, the branches generated by super- and subcritical bifurcations emerge as, 
respectively, stable and unstable ones, and the character of the stability changes when a branch passes a turning 
point. On the other hand, the Vakhitov-Kolokolov criterion^ dE/dk > 0, which in many models with attractive 
nonlinearities is a necessary stability condition [U [3l IH [28] , does not catch the instability of solitons related to the 
symmetry-breaking bifurcations, which is a known fact too [I9l[20]. The stability properties are different for vortices, 
as they may be additionally unstable against azimuthal perturbations [T5HT7] . The analysis of the azimuthal instability 
of vortex solitons in the present model is reported in the next section. 

As seen in Fig. [sj the bifurcation picture for 5 = 1 (without referring to the stability, for the time being) is very 
similar to that for s = 0. The solution branches form a loop, with both the direct and reverse bifurcations being 
subcritical in the interval of < A < 0.0998. At 0.0998 < A < 0.1018, the direct bifurcation is supercritical, while the 
reverse one is still subcritical. The two bifurcations are supercritical, corresponding to the completely convex loop, 
at 0.1018 < A < 0.110 ^ Amax , up to the point where the bifurcation loop ceases to exist. 

The bifurcation loops were also constructed for vortex solitons with spin s = 2. As seen in Fig. [4j the direct and 
reverse bifurcations are subcritical at < A < 0.1001. There is a tiny region (0.1001 < A < 0.1015) in which the 
direct bifurcation is supercritical, while the reverse one stays subcritical. At A 0.1015 the reverse bifurcation also 
switches to the supercritical type, making the loop completely convex. It keeps this shape up to the point of the 
disappearance of the loop, at Amax^^ ^ 0.1102. 



IV. THE LINEAR-STABILITY ANALYSIS 



The stability of stationary solutions ([2| was explored using the standard approach: we take a perturbed solution 
with 



V^(r, 6>, z) = [U{r) + 5U{r, 6>, z)] ex.p{isO) ex.p{ikz), 

(l){r, 6>, z) = [V{r) + SV{r, 6>, z)] exp{isO) exp{ikz), (6) 



where small perturbations are looked for, as usual, in the form of angular eigenmodes, with an integer azimuthal 
perturbation index, n, and the respective instability growth rate 7^ , 



6U = [t/+(r, z) ex.p{inO) + ?7_(r, z) exp(— m6>)] exp(7^z), 

SV = ['l^+(r, z) ex.p{inO) + VL(r, z) exp(— m^)] exp(7nz). (7) 
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FIG. 2. (Color online) The bifurcation diagrams, in the (^totai,©) plane, for fundamental solitons (s = 0), at different values 
of the linear-coupling constant, A. Here a nd in other bifurcation diagrams, stable and unstable branches are shown by solid 
and dotted lines, respectively (see also Fig. 5(a) below, for details of the stability of solution branches displayed in this figure). 
The bifurcation loops produced by the variational approximation (section |v| are shown too, by dashed-dotted red curves. 
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FIG. 3. (Color online) The same as in Fig. [2] but for vortex solitons with s = 1. See also Fig. 6(a) below for details of the 
stability. 



Taking the perturbation in this form leads to a closed system of linearized equations generated by the substitution of 
expressions (|6| and ^ into Eqs. (|3|: 

+(2 - 3U'^)U^U+ + (1 - 2U^)U^U*_ = -XV+; 

-kU- + i-inU- + — ^ + -— ^ 2^U- 

jj^z r or 

+(2 - 3U^)U^U- + (1 - 2U^)U^Ul = -AVI; 

. T. d^V+ ldV+ (s + n)2 

+(2 - 3V^)VV+ + (1 - 2V'^)V'^V1 = -A/7+; 

+(2 - 3V^)V^V- + (1 - 2V^)V^V^ = -XU- . (8) 

These equations are to be solved with the boundary conditions, which demand {U±,V±} r'^^^' at r ^ 0, and the 
exponential decay of the perturbation eigenmodes at r ^ oo. 

There are several available numerical methods for solving such a boundary-value problem and finding the pertur- 
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FIG. 4. (Color online) The same as in Figs. [2] and [s] but for vortex solitons with s = 2. See also Fig. 7(a) below for details of 
the stability. 



bation growth rates, 7^ [Sl [TTJ [26]. We treated Eqs. ([8| as an algebraic eigenvalue problem for 7^, and solved it 
directly, using a finite-difference method. The largest instability-growth rate was identified as the real part of the 
most unstable eigenvalue, max{Re(7n)}. This approach has confirmed the azimuthal stability of the fundamental 
solitons {s = 0) and revealed instability regions for vortices with 5 = 1,2. 

For 5 = 0, the most dangerous perturbation azimuthal index is n = (i.e., as said above, the fundamental solitons 
are not destabilized by azimuthal perturbations). An example of the ensuing curves which display the maximum 
growth rate, in the case corresponding to the bifurcation diagram that features the double bistability at A = 0.05 (see 
Fig. [2|, is shown in Fig. [sj As expected, unstable are backward-going portions of the asymmetric solution branches, 
i.e., in the region between the bifurcation points and turning points of the bifurcation curves. As might be expected 
too, the symmetric solutions are unstable in the entire region between the points of the direct and reverse bifurcations. 

For vortices with s = 1, the growth-rate curves pertaining to n = feature the same behavior as for the fundamental 
solitons. However, additional unstable eigenmodes, for both the symmetric and asymmetric solutions, were found with 
azimuthal indices n = 1, 2, 3. On the other hand, no unstable perturbations were detected for larger n > 3. A typical 
example displaying all the existing unstable eigenvalues for 5 = 1 is presented in Fig. [6| for the same coupling 
constant as in Fig. [sj A = 0.05. In particular, panel (d) in this figure shows the results for the symmetric solutions. 
The eigenmode that remains the last unstable one with the increase of £^totai pertains to n = 2. This instability 
ceases when the energy attains value £^totai ^ 340, which corresponds to propagation constant k ^ 0.2. This value 
is exactly the expected one, according to relation k^^^^"^ ^ 0.15 + A, where k ^ 0.15 is the known stability threshold 
in the single-component model (A = 0), for 5 = 1 [16]. Further, the results for unstable eigenmodes disturbing the 
two inner asymmetric branches (the portions of the branches that commence at the bifurcation point and end at the 
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FIG. 5. (a) The bifurcation diagram, for the fundamental soUtons (s = 0) at A = 0.05, cf. Fig. [5] (b) The corresponding 
maximum growth rate of perturbation eigenmodes for the symmetric and asymmetric solutions. 



turning points) are shown in panelj6[c). At the respective edge points, these curves are linked to their counterparts 
(continuations) in panels |6jd) and|6p)), the latter panel pertaining to the outer asymmetric branches. The conclusion 
is that the inner portions are completely unstable (as well as in the case of 5 = 0, cf. Fig. [5|, while parts of the outer 
branches are stable. In the case shown in Fig. |6] (recall it pertains to A = 0.05), the stability segment of the outer 
branch of the asymmetric solutions with 5 = 1 is 195 < £^totai < 630, which corresponds to 0.165 < k < 0.186. 

With the increase of A, the stable section of the asymmetric branch with s — 1 shrinks, and it disappears at A ~ 0.067 
(which still corresponds to the double-concave shape of the loop with both the direct and reverse bifurcations of the 
subcritical type). At larger values of A, the instability accounted for by the perturbation mode with n = 2 covers the 
entire asymmetric branch. With the further increase of A, the instability regions corresponding to the eigenmodes 
with the other values of n also expand and gradually cover the entire asymmetric branch. In all the cases that we 
have examined, the growth rate corresponding to n = 2 is always the largest one for the upper asymmetric branch. 

In the case of the double vortex, with s = 2, results of the stability analysis are presented in Fig. [7| again for 
A = 0.05. As in the case of 5 = 1, the eigenmode that determines the stability boundaries has n = 2. However, in 
this case the growth rate for perturbations with n = 2 is not necessarily the highest for the outer asymmetric branch, 
and the azimuthal index of the dominant perturbation eigenmode switches from n = 4 to n = 3 and then to n = 2. 
Similar to the behavior of the vortex with s = 1, the stability region of the asymmetric stable solutions diminishes 
with the increase of A. The entire diagram is totally unstable for A > 0.0505, which is slightly larger than A = 0.05 
for which Fig. [7|is displayed. As well as in the case of 5 = 1, the ultimate destabilization occurs when the shape of 
the bifurcation loop is still double-concave. 

Finally, the symmetric vortices with s = 2 are unstable at A < /c < 0.162 + A, where k ~ 0.162 is the stability 
threshold for 5 = 2 in the single-component model [16 . This threshold corresponds to the total energy E ^ 1060 of 
the symmetric vortex. 



V. THE VARIATIONAL ANALYSIS 



The stationary solutions can also be studied analytically by means of the variational approximation (VA), cf. Ref. 
[20] , taking into regard that stationary equations Q can be derived from the Lagrangian, 



-k{U^^V^) 



dU 
dr 



dV 
dr 



-^hu^ + V^) - \{U^ + V^) - 2XUv\ dr. 



(9) 
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FIG. 6. (a) The bifurcation diagram for the vortices with s = 1 at A = 0.05, cf. Fig. [s] The corresponding maximum growth 
rates of perturbation eigenmodes for the inner and outer asymmetric branches and the symmetric one are shown in panels (b), 
(c), and (d), respectively. The labels near the curves indicate the mode's azimuthal number. Note that the plots which appear 
aborted in panels (b) and (c) are actually continuations of each other and of the plots in panel (d). This is in accordance with 
the fact that the outer and inner branches of the asymmetric states are linked at the turning points of the bifurcation diagram, 
and the symmetric and inner asymmetric branches are linked at the bifurcation points. 



To approximate solutions to Eqs. Q (generally, asymmetric ones), the following ansatz was adopted, with common 
width W of both components, but different amplitudes, A and B: 



{[/(r),F(r)}ansatz = R^K ^xp 



(10) 



where 5 = 0, 1,2 is the same spin as above. The energies of the two components of this ansatz, defined according 
to (|4|, are 



The substitution of the ansatz into Lagrangian ([9| and the integration yield the effective Lagrangian: 

2 



(11) 



-Leff = -s\k{A^ ^ B^)W- 



2(l+s) 



{s^l)s\{A'^B')W' 



(^^)' / 44 , d4mt/2(1+2s) 



, (A"" + B'')W'^^'^'^'^ - 

It is convenient to redefine the variational parameters as 



(^6 ^ 56)^2(l+3«) ^ 2s!AABT^2(i+.)^ 



(12) 



a = (W'^/V2j {A + B), p= (wV\/2) {A - B) 
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FIG. 7. The same as in Fig. |6] but for the double vortices (s = 2) at A = 0.05. 



in terms of which effective Lagrangian { 12 ) takes the form of 
2 



-ieff = -(s + l)s!(a2+/^2) 



5!fc(a2 + /?2) + igli(a4+^4^g^2^2) 



4 . 335+2 



(13) 



Values of the variational parameters corresponding to stationary solutions, a, j3 and W, are determined by the 
Euler-Lagrange equations, dLefi/d{W^) = dLes/d{a'^) = dLes/d{l3'^) = 0, i.e., 



(2^)! 



(3s)! 



4 • 33^+2s! 



r2 , ,^ , (2s)! (c,2 + 3;j2)^2 



-fcl^^ - (s + 1) + 



22s+2gl 



(3s)! 



4 ■ 33^+is! 
- _ _^ 1) _^ 

4.33.+I5!' 



(a^ + 10a2/?2 + 5;9^)W^2 ^ ^^^2 ^ g 
(/32 + 3a2)H^2 



22«+2s! 

+ lOa^/?^ + 5a'*)M^2 - = 0. 



(14) 



(15) 



(16) 
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The variational solutions are obtained by numerically solving the system of equations, (14)-(16) for a, /3 and 



for given 5, A and k. In this way, several VA-predicted bifurcation loops were constructed, for different values of A 
and for 5 = 0, 1 and 2, as shown above in Figs. [2j [3] and [4] by dashed-dotted curves, alongside the numerically found 
loops. The figures show that the VA quite accurately predicts the transformation of the bifurcation loop from the 
concave shape to the convex one with the increase of A. An adequate indication of the accuracy of the VA is given 
by comparing critical values of A at which the symmetry-breaking bifurcations disappear, along with the bifurcation 
loops. For that purpose, we set /3 = in Eqs. (14)-(16) and subtract the second equation from the third, which yields 



2 33^+15! 



(3s)! I 22^+2s! 



± 



228+2^1 



^335+1^! 



A 



(17) 



Within the framework of of the VA, the asymmetric solutions exist under the condition that expression (17) yields 
real values: 



A < A.s 



33.+i^2g)!)2 
24^+5(35)!s! ' 



(18) 



For 5 = 0, 1 and 2, Eq. (18) predicts critical values Aq = 0.09375, Ai = 0.10547, and A2 = 0.10679. The comparison 



with their numerically found counterparts shows that the difference is 2.8%, 4.3% and 3.9%, respectively. 



To calculate coordinates of the VA-predicted bifurcation points, we substitute expression (17) into Eq. (16) with 



/3 = 0. Both the variational and the numerically generated plots for values of E and k at the bifurcation points are 
shown, versus the coupling constant. A, in Fig. |8] It is seen that the variational and numerical results are always in 
good agreement for the direct bifurcation. On the other hand, the approximation for the reverse bifurcation becomes 
inaccurate for very small values of A. This difference is explained by the fact that, near the reverse bifurcation, the 
actual profiles of the soliton components become increasingly rectangular-like, i.e., different from the shape assumed 
by ansatz (10). 



VI. DEVELOPMENT OF THE INSTABILITY OF VORTEX RINGS 

To explore results of the instability development, direct simulations of Eqs. ([T]) were performed by means of the 
standard pseudospectral split-step method, for initial conditions corresponding to the stationary solutions presented 
in section [nil Perturbations were not explicitly added to unstable solitons, the instability being initiated by truncation 
errors of the numerical code. First, we demonstrate the splitting of vortex solitons which are unstable against azimuthal 
perturbations. For the single-component model, a similar numerical analysis was reported in Ref. [T6|, where it was 
concluded that, generally, the azimuthal index, n, of the most unstable eigenmode determines the number of fragments 
produced by the splitting. 

Figure [9] displays the numerically simulated evolution of the asymmetric vortex ring with 5 = 1, A = .05, k = 0.155 
and £^totai ^ 150, for which the single unstable perturbation eigenmode has n = 2, as per Fig. |6(b)[ In this case, 
the breakup of the vortex becomes conspicuous at Zspia ^ 900, giving rise to two fragments, in accordance with the 
linear-stability analysis. 



Similar results for the double asymmetric vortices {s = 2) and A = 0.05 are presented in Figs. 10 -12 The initial 
states were chosen so as to have, in each case, the largest growth rate at a different value of n. To this end, we took 
k = 0.112, 0.1505, 0.174, which correspond to the vortices with £^totai ^ 165, 250, 500, the corresponding largest 
instability growth rates being 7n=4 ^ 0.098, 7n=3 ^ 0.047, 7n=2 ^ 0.009, respectively. As expected, the numbers of 
fragments generated by the breakup are consistent with these values of n. We stress that, in all the cases presented 
here, there is a well-pronounced dominant eigenmode. As mentioned in Ref. [16 , when the parameters are taken 
close to borders between regions dominated by unstable eigenmodes with different azimuthal indices n, it is difficult 
to predict which one will determine the outcome of the splitting. 

Values of the propagation distance needed for the splitting to commence are also in agreement with the predictions 
based on the growth rates of the linear instability. In particular, a large growth rate was found for A: = 0.112 (£^totai ^ 
165) and, accordingly, in that case the splitting starts very early, at Zspia ~ 170. For k = 0.1505 (^totai ^ 250), the 
breakup starts later, at ^spiit ^ 680, and when the growth rate is small - for instance, at k = 0.174 (£^totai ^ 500) - 
the splitting sets in after a very long evolution, at 2;spiit ^ 3250. In all the cases that we have examined, the fragments 
maintain the asymmetry of the original unstable vortex rings. 

Note that unstable asymmetric solutions could transform into stable symmetric ones (and vise versa) if the only 
unstable perturbation eigenmode were the one with n = 0. In fact, this happens solely for 5 = (see Fig. [sj for 
example). In all the cases that we have examined, the azimuthal instability of the vortices with 5 = 1 and 2 destabilizes 
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XXX 

(a) (b) (c) 




XXX 

(d) (e) (f) 



FIG. 8. The comparison of the numericahy generated (circles connected by continues hnes) and variationally predicted (dashed- 
dotted hnes) points of the direct and reverse bifurcations, (a)-(c) The propagation constant, k, at which the bifurcations occur, 
as a function of A, for s = 0, 1 and 2. (d)-(f) The same for the total energy, ^totai- 

and destroys the solutions, before they could be reshaped into stable symmetric or asymmetric structures with the 
same s. 



VII. INTERACTIONS BETWEEN SOLITONS 



Direct simulations were also used to study interactions between two initially quiescent solitons separated by a 
relatively small distance. In the 2D single-component model, a similar investigation was reported in Ref. [14j, for 
vortices with 5 = 1. Collisions between asymmetric solitons in the ID dual-core model were studied earlier in Ref. [29] . 
Here we focus on pairs of initially quiescent solitons (rather than moving ones), as the interaction effects are strongest 
in such a case. 

Our analysis was performed for several combinations of stable asymmetric solutions, for all the three values of the 
spin considered here, s = 0, 1,2. First, we examined the interaction between identical solitons separated by distance 
Ax, namely: {U,V}.^.^.^^{x,y) = ^Istationary - ^) + {^^ ^Istationary + ^cxt, wc cousidcrcd 

pairs of cross-identical asymmetric solitons (one being a specular counterpart of the other): {U,V}^^^^^^^{x^y) = 
{t/, "^Istationary ~ v) + ^istationary + v) ' ^dive also performed the simulatious for the soliton pairs 
with the phase shift of = tt. 

The results are presented in Figs. p!3||T6| In accordance with the known principle [30 , we observed that identical 
solitons with even values of the spin, s = or 5 = 2, attract each other in the in-phase configuration, while the 
vortices with 5 = 1 exhibit the attrac t ion when t hey are out-of-p hase^ with A6 = tt. The attraction results in inelastic 
collisions, as seen in Figs. [13(7 " 



(b)| [T6(a)jj(b)| and [T5(a)jj(b)| If both solitons are fundamental ones, the inelastic 



interaction ends up with their merger into a single pulse (which is not necessarily another stationary fundamental 
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FIG. 9. Gray-scale plots illustrating the splitting of an un- 
stable vortex solution, at s = 1 and k = 0.155 [Etotai ~ 150). FIG. 10. The same as in Fig.[9j but for s = 2 and k — 0.174 
Values of the propagation distance are labeled above the left (Etotai ~ 500). 
frames. 



soliton). If vortex rings are involved into the collision, the eventual result is destruction of the ring(s), and formation 
of a disordered pattern. 

On the other hand, repulsion is observed, also in agreement with the predictions of Ref. [30^, between out-of-phase 
solitons (A^ = tt) with even values of the spin, 5 = and 2 (not shown here in detail), and between in-phase vortices 
with 5 = 1. In the case of the repulsion, the interactions produce, as a matter of fact, no visible eff ect, leading only 
to a small increase of the separation between the two solitons, as shown, for s = 1, in Fig |14(a)] - |(b)| 

Additional simulations were carried out for soliton pairs with different spins. Figure [T7| shows a typical example 
for an in-phase vortex pair with s = 1 and s = 2. In this case, the phase varies smoothly between the centers of 
the two vortices, similar to the case of the in-phase pairs with 5 = or s = 2, and the out-of-phase one with 5 = 1. 
Accordingly, the simulations demonstrate a strong attractive interaction in this case. On the other hand, for the 
same pair with 5 = 1, 5 = 2 and = tt, the change of the phase between the vortices is very steep, similar to the 
situations when the solitons with equal spins repel each other, and, accordingly, the vortices under consideration repel 






each other too, which results in slow separation between the vortices (not shown here). 

We have also studied the interaction between the vortex with 5 = 1 and the fundamental soliton (s = 0). In this 



case, an inelastic collision is observed in Fig. [18] for A6> = tt, while for A6> = the interaction is repulsive, producing 
no conspicuous effect (not shown here). These outcomes may be explained by the same character of the phase pattern 
between the solitons as above, i.e., smooth in the configuration leading to the attraction and strong interaction, and 
steep in the opposite case of the repulsion. 



Interactions between stable symmetric solutions have also been examined, demonstrating results identical to those 
originally reported in Ref. [14] (not shown here). In particular, in this case the collisions maintain the symmetry 
between the two components, and do not lead to the appearance of any asymmetric final states. 
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FIG. 13. The interaction between fundamental (s = 0) in-phase (A^ = 0) asymmetric solitons, with A = 0.05, k = 0.16, 
and initial separation Ax = 30. In panel (a) the solitons are identical, while in (b) they are cross-symmetric. The evolution 
distances, z, are marked above each frame. 



VIII. CONCLUSION 



We have introduced a model of a 2D dual-core waveguide with the CQ (cubic-quintic) nonlinearity inside each core 
and the linear coupling between them. Families of fundamental {s = 0) and vortical solitons, with spins 5 = 1 and 
2, have been constructed, and their stability has been investigated. The model may be realized as a dual-core planar 
optical waveguide, or as a set of two tunnel-coupled parallel pancake- shaped traps for BEG. In the former case, the 
solitons are "planar light bullets" . In particular, vortex solitons may be interpreted as a new species of the "bullets" , 
viz., spatiotemporal vortices. 

The main objective of the work was to study symmetry-breaking bifurcations of the 2D solitons, both fundamental 
and vortical ones. If the inter-core coupling constant. A, is not too large, the bifurcation diagrams for the solitons 
of all the types form closed loops, which connect points of direct and reverse bifurcations. The stability of all 
the solutions was investigated via the calculation of the corresponding eigenvalues for infinitely small perturbations 
around the solitons. In particular, it was found that, at sufficiently small values of A, the loop for the fundamental 
solitons {s = 0) includes two bistability regions, which may be of interest for potential applications, such as all-optical 
switching. This is also a new feature of the loop in comparison with its earlier studied counterpart in the ID version of 
the model, which could feature a single bistability domain. At larger values of A, both bistable regions vanish and the 
loop's shape becomes plainly convex. With the further increase of A, the loop shrinks to zilch and disappears, leaving 
only stable symmetric solitons. The vortical solitons may be easily destabilized by azimuthal perturbations, but they 
also have stability regions, as long as the corresponding bifurcation loop keeps its double-concave shape. We have 
also developed a quasi- analytical approach to the description of the bifurcation diagrams, based on the variational 
approximation, which produces reasonably accurate predictions, in comparison with the numerical results. 

In direct simulations, we have demonstrated the splitting of azimuthally unstable asymmetric solitons into sets of 
fragments. The number of the fragments usually corresponds to the azimuthal index of the most unstable eigenmode 
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FIG. 14. The interaction between asymmetric in-phase vortices with s — 1, A = 0.05, k — 0.18, and Ax = 40. Notice the 
repulsion between the in-phase vortex sohtons in this case. 

of small perturbations. We have also studied interactions between initially quiescent solitons, and confirmed the 
earlier prediction [30 , which states that the usual attractive/repulsive sign of the interaction between the in-phase/7r- 
out-of-phase solitons with even values of the spin (5 = or 2), is reversed for the odd spin (5 = 1). In the case of the 
attraction, the vortex solitons merge into disordered patterns, losing the initial topological structure. 

This work may be naturally extended in other directions. In particular, it may be interesting to study symmetry- 
breaking effects in two-component 2D solitons and vortices in the system of NLSEs coupled by both linear and 
nonlinear terms, which may describe the co-propagation of two polarizations of light in a single nonlinear waveguide 
[1]. Moreover, the latter version of the model is meaningful in the 3D geometry too. Another relevant generalization 
still pertains to the dual-core waveguide, with the linear coupling between the two waves, while the competing nonlinear 
terms are quadratic and cubic, rather than cubic and quintic, cf. Ref. [31 and references therein. In that case, one 
may also expect bifurcation loops accounting for the breaking and restoration of the symmetry of two-component 
solitons. 
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FIG. 15. The interaction between asymmetric vortices with s = 1, = tt, A = 0.05, k = 0.18, and Ax = 40. Notice the 
attraction between the 7r-out-of-phase vortex sohtons in this case. 
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FIG. 16. The interaction between asymmetric in-phase vortices with s = 2, A = 0.05, k = 0.184, and Ax = 60. 
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FIG. 17. The interaction between vortices with s = 1 and 
s = 2, for AO = 0, Ax = 50, and A = 0.05, k = 0.184. 
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FIG. 18. The interaction between a fundamental sohton 
(s = 0) and a vortex with s = 1, for Ax = 34, = tt and 
A = 0.05, k = 0.17. 



